Elena Spielmann, Patrick Jones, Leanne Chook, and Maeve Grady
Introduction
Obesity has reached epidemic proportions in the United States, and affects people of all ages, socioeconomic backgrounds and ethnicities. This imposes substantial economic burdens in the form of productivity losses, strains to the healthcare system and costs.
The prevalence of obesity in the United States was approximately 41.9% in 2017 to 2020, of that, the prevalence of severe obesity also increased to 9.2%. This increasing trend is worrying because obesity can lead to several other health conditions such as heart disease, strokes, type 2 diabetes, and certain types of cancers, all of which are the leading causes of preventable and premature deaths in the country. For these reasons, determining predictors of obesity is a public health priority.
Furthermore, this increasing trend of obesity also has other implications on national security, especially as more individuals become ineligible for the military, therefore, significantly reducing military recruitment numbers.
Motivating question:
What health, social, and economic indicators are most important in predicting obesity?
The main motivating question for this project was to examine obesity rates across the United States and find predictors that can help policymakers make effective decisions to help reduce obesity prevalence. Each of the selected data sets brings in different predictors to the model, covering health, social and economic aspects of obesity. This is relevant to the public policy sphere because the marrying of these different aspects are key to meaningful policy change and effective spending of government funds. Determining important predictors can help policymakers make effective decisions to help reduce obesity prevalence, which is a public health priority.
Outline of project:
We will examine and clean up the available information from datasets downloaded from PLACES and the Food Access Repository Atlas (FARA). We will then conduct exploratory data analysis (EDA) and geospatial analysis to guide our decision for the final models, before creating several predictive models to find the best model with highest predictive power. We will also determine what the most important variables for predicting obesity are for each model.
This endeavor used both the CDC’s PLACES dataset and the USDA’s FARA dataset. The team attempted to read each dataset in using API query’s (attempts below) but ran into limitations in the APIs and in technical expertise. In the case of the FARA dataset, the API available is hosted through ESRI as an ArcGIS REST service. This means that the query must be read in through an st_read()command, which made it very difficult to troubleshoot query errors. There was limited documentation available for the FARA API.
After joining these two datasets, we also read in census tract shapefiles using the census API and library(tidycensus).
Code
# | eval: false# accessing FARA data through ArcGIS REST API## to write this code I consulted this blog post: https://community.esri.com/t5/gis-blog/accessing-arcgis-rest-services-using-r/ba-p/898451url <-parse_url("https://gis.ers.usda.gov/arcgis/rest/services")url$path <-paste(url$path, "foodaccess2019/MapServer/0/query", sep ="/")url$query <-list(returnGeometry ="true",f ="geojson",outFields ="*")request <-build_url(url)foodaccess <-st_read(request)
Reading layer `GeoJSONSeq' from data source
`https://gis.ers.usda.gov/arcgis/rest/services/foodaccess2019/MapServer/0/query?returnGeometry=true&f=geojson&outFields=%2A'
using driver `GeoJSONSeq'
Simple feature collection with 1 feature and 0 fields (with 1 geometry empty)
Geometry type: GEOMETRYCOLLECTION
Dimension: XY
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: WGS 84
Code
## using this api to try to figure out how the query should work:## https://gis.ers.usda.gov/arcgis/rest/services/foodaccess2019/MapServer/0/query?where=&text=&objectIds=&time=&geometry=&geometryType=esriGeometryPolygon&inSR=&spatialRel=esriSpatialRelIntersects&distance=&units=esriSRUnit_Foot&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&havingClause=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false&resultOffset=&resultRecordCount=&returnExtentOnly=false&datumTransformation=¶meterValues=&rangeValues=&quantizationParameters=&featureEncoding=esriDefault&f=html
Code
#install.packages("RSocrata")library(RSocrata) ## the Rsocrata package has to be loaded in after httr package for the above code to work ## reading in api credentialsload_dot_env(here(".env"))app_token <-Sys.getenv("PLACES_app_token")user_name <-Sys.getenv("PLACES_username")password <-Sys.getenv("PLACES_password")# PLACES datasetplaces <-read.socrata("https://chronicdata.cdc.gov/resource/swc5-untb.json",app_token =paste(app_token),email =paste(user_name),password =paste(password))
Combining the Datasets
The federal government provides a wide variety of publicly available health outcomes data. Our group explored a variety of different sources, but ultimately settled on two datasets: the USDA’s The Food Access Research Atlas and the CDC’s PLACES database.
The Food Access Research Atlas (FARA) provides a variety of food access measures for low income and low access census tracts. Food access measures use income, transportation, and distance from grocery stores and other food sellers to determine how accessible food is to residents of a given census tract. This dataset includes distinct measures of access for urban and rural populations.
PLACES provides a variety of health outcome measures including obesity, arthritis, and diabetes, as well as health behaviors such as sleep and smoking. The data in places is largely drawn from the CDC’s Behavioral Risk Factor Surveillance system survey and the National Survey of Children’s health. Data from PLACES is organized at the census tract level.
Our team downloaded both data files directly from the CDC/USDA websites. Because we worked with two separate datasets, our team had to merge these two datasets together before conducting exploratory data analysis and modeling. The USDA organizes the FARA dataset in the tidy format, the dataset needed minimal cleaning and manipulation prior to merging. On the other hand, the CDC does not organize the PLACES dataset in the tidy format. In its original format, PLACES lists data in the long format, with each variable individually listed for each census tract. Our team pivoted this data to the wide format, giving each variable its own column. We also dropped a number of variables that were not relevant to our analysis, such as data source, footnotes, and upper and lower limit estimates. After converting the PLACES dataset, our team successfully merged both datasets together.
Our team wanted to perform geospatial exploratory data analysis in addition to standard EDA. To accomplish this, we pulled census tract shape files from the CDC using their publicly available API. We merged this shapefile to our combined dataset to enable geospatial EDA, we set this combined dataset’s CRS to 4326.
Code
# load the FARA dataset fara <-read_excel("fara_2019.xlsx")# clean the variable names fara <-clean_names(fara)# load the PLACES datasetplaces <-st_read("PLACES_ Local Data for Better Health, Census Tract Data 2022 release.geojson") %>%st_transform(value =4326)
Reading layer `PLACES_ Local Data for Better Health, Census Tract Data 2022 release' from data source `/Users/leannechook/Desktop/school_documents/gradschool/spring_2023/data_science/assignments/ppol670_finalproject/PLACES_ Local Data for Better Health, Census Tract Data 2022 release.geojson'
using driver `GeoJSON'
Simple feature collection with 2161543 features and 22 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -173.0177 ymin: 19.25282 xmax: -67.00993 ymax: 71.23014
Geodetic CRS: WGS 84
Code
places <- places %>%rename(census_tract = locationid)# creating tidy datasetplaces_tidy <- places# Filtering observations to 2020 onlyplaces_tidy <- places %>%filter(year ==2020)# There are a number of columns that are not necessary, they need to be dropped as wellplaces_tidy <- places_tidy %>% dplyr::select(-(c(statedesc, datasource, category, measure, data_value_unit, data_value_type, low_confidence_limit, high_confidence_limit, categoryid, datavaluetypeid, short_question_text, countyfips, locationname, data_value_footnote, data_value_footnote_symbol, )))# Finally, we can pivot this dataset to wide formatplaces_tidy <- places_tidy %>%pivot_wider(names_from ="measureid",values_from ="data_value")# merge the data frames by census tract# swapped places and fara so that the combined dataset only has census tracts that are in faracombined <-right_join(places_tidy, fara,by ="census_tract")# Read in the shapefile - opting for tidy census heredotenv::load_dot_env()api_key <-Sys.getenv("CENSUS_API_KEY")options(tigris_use_cache =TRUE)allstates <-c(state.abb)allstates
census <-get_acs(geography ="tract", variables ="B19013_001",year =2019,state = allstates, geometry =TRUE) census <- census %>%rename(census_tract = GEOID) %>%st_transform(crs =st_crs(4326))# binding census shapefile to this datasetcombined_nonsf <- combined %>%st_transform(crs =st_crs(4326))combined_census <-st_join(census, combined_nonsf, join = st_intersects)
Exploratory Data Analysis
Before building models, our team conducted exploratory data analysis in order to identify variables of interest and their relationships with our outcome variable (obesity). An initial glimpse at our dataset revealed that a large number of numeric variables were coded as character variables, making them unsuitable for analysis and modeling. We identified and converted these variables to numeric format, while leaving character variables untouched. We recoded a number of binary variables as numeric rather than factors, as factors interfered with our models.
Code
# A ton of our numeric variables are stored as character variables, they need to be converted to numeric variablesto_convert <-colnames(select_if(combined, is.character))to_convert <- to_convert[!to_convert %in%c("census_tract", "countyname", "year", "stateabbr", "geometry")]to_convert <-unlist(to_convert)combined <- combined %>%mutate_at(to_convert, as.numeric)
Comparing Obesity Prevalence with Variables of Interest
We continued exploring our data by creating a number of graphs of the relationships between variables of interest and obesity. The variables of interest that we focused on are: 1) Poverty Rate; 2) Median Family Income; 3) Prevalence of Poor Mental Health; and 4) Frequency of Health Checkups.
To make this process easier, we wrote custom functions to make multiple plots.
Code
#' eda_point documentation#'#' @param x An independent numeric variable such as poverty rate#' @param xlab A human readable title for the x-axis#' @return A geom_point scatter plot of x vs y by fill#' @export Does not export#'#' @examples#' eda_point <-function(x, xlab) { plot <-ggplot(data = combined, aes(x = {{ x }},y = OBESITY)) +geom_point(aes(fill =as.factor(urban)),color ="white", shape =21) +labs(x = xlab, y ="Obesity Prevalence") +scale_fill_discrete(name =NULL,breaks =c(0, 1),labels =c("Rural", "Urban")) +theme_minimal()return(plot)}eda_point(poverty_rate, "Poverty Rate")
Code
eda_point(median_family_income, "Median Family Income")
Code
eda_point(MHLTH, "Prevalence of Poor Mental Health")
Code
eda_point(CHECKUP, "Frequency of Health Checkups")
These plots reveal a number of insights about our data. While the average obesity rate is higher in rural areas, urban areas have higher variability of data. Unsurprisingly, as the poverty rate increases, obesity rates increase as well. The next graph reflects this as well, showing that obesity declines as median family income increases. As expected, poor mental health also has a positive correlation with obesity prevalence. Regular health checkups, however, do not have a clear relationship with obesity. Furthermore, there is a much wider dispersion of rural census tracts in this plot as well.
Simple feature collection with 6 features and 4 fields
Geometry type: MULTIPOINT
Dimension: XY
Bounding box: xmin: -173.0177 ymin: 30.2377 xmax: -84.99064 ymax: 71.23014
Geodetic CRS: WGS 84
# A tibble: 6 × 5
stateabbr mean_obesity median_income sd_food_access geometry
<chr> <dbl> <dbl> <dbl> <MULTIPOINT [°]>
1 AK 32.1 84524 0.501 ((-156.7348 71.23014), (-…
2 AL 40.3 54489 0.495 ((-85.58971 34.58632), (-…
3 AR 38.0 53672 0.500 ((-90.86605 34.31843), (-…
4 AZ 31.0 62880 0.497 ((-113.382 36.52239), (-1…
5 CA 29.1 78986 0.430 ((-120.5514 41.50297), (-…
6 CO 25.3 80828. 0.489 ((-102.6192 38.09705), (-…
Code
#' eda_bar#'#' @param y A variable of interest to be plotted by state #' @param ylab A human readable title for the y axis#'#' @return#' @export#'#' @examples#' eda_bar <-function(y, ylab){ plot <-ggplot(summary_stats, aes(x = stateabbr, y = {{ y }})) +geom_bar(stat ="identity", fill ="dodgerblue") +labs(x ="State", y = ylab) +theme(axis.text.x =element_text(angle =90, hjust =1)) +theme_minimal()return(plot)}eda_bar(mean_obesity, "Mean Obesity")
Code
eda_bar(median_income, "Median Income")
Code
eda_bar(sd_food_access, "Standard Deviation of Food Access")
Code
# boxplot to compare the distribution of obesity prevalence between urban and rural areasggplot(combined, aes(x = urban, y = OBESITY)) +geom_boxplot()
Spread of obesity prevalance across the United States
Code
combined_census <- combined_census %>%clean_names() %>%filter()breaks <-c(0, 10, 20, 30, 40, 50, 60, 70)# binning obesity rates and creating text to work withcombined_census <- combined_census %>%mutate(obesity_num =as.numeric(obesity),obesity_bin =cut(obesity_num, breaks =10),maptext =paste0("State: ", as.character(stateabbr), "\n","Tract: ", as.character(census_tract_x), "\n","Obesity rate: ", as.character(obesity), "\n","Rate of current smokers: ", as.character(csmoking), "\n","Rate of teeth lost: ", as.character(teethlost), "\n","Median family income: ", as.character(median_family_income)))# plot the graph obesity_tracts <-ggplot(combined_census) +geom_sf(aes(fill = obesity_bin),lwd =0,color =NA)+scale_fill_brewer(palette ="RdYlGn", direction =-1)+labs(title ="Adult obesity rate by census tract in the US",fill ="Prevalence of\nObesity") +coord_sf(xlim =c(-160, -64),ylim =c(18, 72), crs =4326) +theme_void()obesity_tracts
Code
#this codeblock is an interactive map of obesity by census tract that includes some other information that can pop up when hovering over each tract. While the code runs, it takes a very long time to run so we have chosen not to evaluate it here. interactivemapdata <- combined_census %>%select(name, state, census_tract_x, obesity_num, median_family_income, csmoking, teethlost, maptext) interactivemapdata %>%mapview(zcol ="obesity_num",map.types ="Esri.NatGeoWorldMap",col.regions =brewer.pal(8, "RdYlGn"),color =NA)obesity_tracts_int <-ggplot() +geom_sf_interactive(data = combined_census, aes(fill = obesity_bin,data_id = maptext,tooltip = maptext),lwd =0,color =NA) +scale_fill_brewer(palette ="RdYlGn",direction =-1, aesthetics ="fill", guide ="none")+labs(title ="Adult obesity rate by census tract in the US")+coord_sf(xlim=c(-160, -64),ylim=c(-30,72), crs =4326) +theme_void()girafe(ggobj = obesity_tracts_int) %>%girafe_options(opts_hover(css ="fill:blue;"), opts_sizing(rescale =TRUE))
(Supervised) Machine Learning
Setting up the testing environment
Code
# separate geometry column into longitude and latitude, and reduce the number of variables to prepare it for recipe combinedsmall <- combined %>%select(census_tract, COPD, OBESITY, STROKE, DEPRESSION, LPA, CASTHMA, TEETHLOST, ARTHRITIS, DIABETES, BINGE, SLEEP, ACCESS2, PHLTH, DENTAL, MHLTH, CANCER, CHD, GHLTH, CHECKUP, CSMOKING, CERVICAL, KIDNEY, COLON_SCREEN, COREW, urban, median_family_income, la1and10, lalowi1share, lakids1share, laseniors1share, lahunv1share)combinedsep <- combinedsmall %>%mutate(longitude =st_coordinates(geometry)[, "X"],latitude =st_coordinates(geometry)[, "Y"]) # to remove the geometry column from dataframe combinedsep <-st_drop_geometry(combinedsep)# convert all variables to numericcombinedsep$urban <-as.numeric(as.character(combinedsep$urban))# set seedset.seed(20230507)# split the data into training and testing sets obesity_split <-initial_split(data = combinedsep, prop =0.8)obesity_train <-training(x = obesity_split)obesity_test <-testing(x = obesity_split)# set up v-fold cross validation folds <-vfold_cv(data = obesity_train, v =5, repeats =1)# create a recipe obesity_rec <-recipe(OBESITY ~., data = obesity_train) %>%step_other(census_tract) %>%step_naomit(all_predictors()) %>%step_center(all_numeric_predictors()) %>%step_scale(all_numeric_predictors()) %>%step_dummy(census_tract)
Models
Linear Regression
This is a simple and widely-used statistical model for predicting a continuous outcome based on one or more predictor variables. The goal is to find a linear relationship between the outcome and predictors, which can be used to make predictions.
Code
# create a linear regression modellm_mod <-linear_reg() %>%set_engine("lm")# create a workflow lm_wf <-workflow() %>%add_recipe(obesity_rec) %>%add_model(lm_mod) # fit the model lm_cv <- lm_wf %>%fit_resamples(resamples = folds)# select the best model based on the "rmse" metriclm_best <- lm_cv %>%select_best("rmse")# finalize workflowlm_final <-finalize_workflow( lm_wf,parameters = lm_best)# fit to the training data and extract coefficientslm_coefs <- lm_final %>%fit(data = obesity_train) %>%extract_fit_parsnip() %>%vi(lambda = lasso_best$penalty)
LASSO Model
This is a linear regression model with a regularization term added to the objective function, which helps to prevent overfitting and improve generalization performance. The LASSO penalty shrinks some of the regression coefficients towards zero, effectively selecting a subset of the most important predictors for the outcome.
Code
# create a tuning grid for lasso regularization, varying the regularization penaltylasso_grid <-grid_regular(penalty(), levels =10)# create a linear_regression model to tune the penalty parameterlasso_mod <-linear_reg(penalty =tune(), mixture =1) %>%set_engine("glmnet")# create a workflow using your updated linear regression model lasso_wf <-workflow() %>%add_recipe(obesity_rec) %>%add_model(lasso_mod) # perform hyperparameter tuning using the lasso_grid and cross validation foldslasso_cv <- lasso_wf %>%tune_grid(resamples = folds,grid = lasso_grid )# select the best model based on the "rmse" metriclasso_best <- lasso_cv %>%select_best(metric ="rmse")# finalize workflow and finding the best modellasso_final <-finalize_workflow( lasso_wf,parameters = lasso_best)# fit to the training data and extract coefficientslasso_coefs <- lasso_final %>%fit(data = obesity_train) %>%extract_fit_parsnip() %>%vi(lambda = lasso_best$penalty)
Ridge Model
This is also a linear regression model with a regularization term, but instead of selecting a subset of predictors, it shrinks all of the regression coefficients towards zero. This can help to reduce the impact of multicollinearity among the predictors.
Code
# create a tuning grid for ridge regularization, varying the regularization penaltyridge_grid <-grid_regular(penalty(), levels =10)# create a linear_regression model to tune the penalty parameterridge_mod <-linear_reg(penalty =tune(), mixture =0) %>%set_engine("glmnet")# create a ridge workflow using your updated linear regression model ridge_wf <-workflow() %>%add_recipe(obesity_rec) %>%add_model(ridge_mod)# perform hyperparameter tuning using the on ridge hyperparameter grid and use cross_validation folds ridge_cv <- ridge_wf %>%tune_grid(resamples = folds,grid = ridge_grid )# select the best model based on the "rmse" metricridge_best <- ridge_cv %>%select_best(metric ="rmse")# finalize workflowridge_final <-finalize_workflow( ridge_wf,parameters = ridge_best)# fit the final ridge model to the full training data and extract coefficientsridge_coefs <- ridge_final %>%fit(data = obesity_train) %>%extract_fit_parsnip() %>%vi(lambda = ridge_best$penalty)
Elastic Net
This is a combination of LASSO and ridge regression, with both penalties included in the objective function. This can help to balance the benefits of variable selection and coefficient shrinkage, and is useful when the data contains many correlated predictors.
Code
# create a tuning grid for elastic net regularization, varying the regularization penaltyelastic_net_grid <-grid_regular(penalty(), levels =10)# create a linear_regression model to tune the penalty parameterelastic_net_mod <-linear_reg(penalty =tune(), mixture =0.5) %>%set_engine("glmnet")# create an elastic net workflow using your updated linear regression modelelastic_net_wf <-workflow() %>%add_recipe(obesity_rec) %>%add_model(elastic_net_mod)# perform hyperparameter tuning using the on your elastic net hyperparameter grid and use cross_validation folds elastic_net_cv <- elastic_net_wf %>%tune_grid(resamples = folds,grid = elastic_net_grid)# select the best model based on the "rmse" metricelastic_net_best <- elastic_net_cv %>%select_best(metric ="rmse")# finalize workflowelastic_net_final <-finalize_workflow( elastic_net_wf,parameters = elastic_net_best)# fit the final elastic net model to the full training data and extract coefficientselastic_net_coefs <- elastic_net_final %>%fit(data = obesity_train) %>%extract_fit_parsnip() %>%vi(lambda = elastic_net_best$penalty)
Comparing the Different Models
Code
# the models are compared for prediction accuracybind_rows(`lm`=show_best(lm_cv, metric ="rmse", n =1),`LASSO`=show_best(lasso_cv, metric ="rmse", n =1),`ridge`=show_best(ridge_cv, metric ="rmse",n =1),`enet`=show_best(elastic_net_cv, metric ="rmse", n =1),.id ="model")
# A tibble: 4 × 8
model .metric .estimator mean n std_err .config penalty
<chr> <chr> <chr> <dbl> <int> <dbl> <chr> <dbl>
1 lm rmse standard 2.55 5 0.00754 Preprocessor1_Model1 NA
2 LASSO rmse standard 2.55 5 0.00747 Preprocessor1_Model01 1e-10
3 ridge rmse standard 2.70 5 0.00735 Preprocessor1_Model01 1e-10
4 enet rmse standard 2.55 5 0.00749 Preprocessor1_Model01 1e-10
Code
all_coefs <-bind_rows(`lm`= lm_coefs,`LASSO`= lasso_coefs,`ridge`= ridge_coefs,`enet`= elastic_net_coefs,.id ="model") all_coefs %>%group_by(model) %>%slice_max(Importance, n =10) %>%ggplot(aes(Importance, Variable, fill = model)) +geom_col(position ="dodge")
Code
all_coefs %>%filter(model !="lm") %>%group_by(model) %>%slice_max(Importance, n =10) %>%ggplot(aes(Importance, Variable, fill = model)) +geom_col(position ="dodge")
Code
# compare the regularized coefficients to the lm coefficients for all three modelsplot1 <-left_join(rename(lm_coefs, lm = Importance),rename(lasso_coefs, LASSO = Importance),by ="Variable") %>%ggplot(aes(lm, LASSO)) +geom_point(alpha =0.3)plot2 <-left_join(rename(lm_coefs, lm = Importance),rename(ridge_coefs, ridge = Importance),by ="Variable") %>%ggplot(aes(lm, ridge)) +geom_point(alpha =0.3)plot3 <-left_join(rename(lm_coefs, lm = Importance),rename(elastic_net_coefs, enet = Importance),by ="Variable") %>%ggplot(aes(lm, enet)) +geom_point(alpha =0.3)plot1 + plot2 + plot3
Decision Tree Model
This is a tree model. The first model is not tuned and the second model is tuned using the same cross-folds as the above models, but does not show an improvement in rmse compared to the non-tuned model. This model was interesting to examine to understand what factors the model considers as driving obesity rates more than others. As you can see in the variable importance plot below, current smoking rates, lost teeth rates, “general health”, low rates of no liesure-time physical activity and physical health as relatively important variables in predicting obesity rates. Median family income was also in the top 10 of important variables.
Code
# create a modeltree_mod <-decision_tree() %>%set_engine(engine ="rpart") %>%set_mode(mode ="regression")# create a workflowtree_wf <-workflow() %>%add_recipe(obesity_rec) %>%add_model(tree_mod)# fit model to the training set tree_fit <- tree_wf %>%fit(data = obesity_train)# create a tree rpart.plot::rpart.plot(x = tree_fit$fit$fit$fit, roundint =FALSE)
In order to pick the best model, we consider the RMSE as the metric used for evaluation, the mean value of the metric across the folds, and the standard error of the metric. Generally, we want to choose the model that has the lowest value of the metric, as this indicates better performance on the test data.
The RMSE values for lm, LASSO, and enet are all the same at 2.55. However, the ridge model has a slightly higher RMSE of 2.69 and the decision tree model has an rmse of 3.96. This suggests that the lm, LASSO, and enet models may be better choices than the ridge model.
Additionally, we consider the standard error of the metric. The standard error gives an estimate of the precision of the mean estimate, with a lower value indicating more precision. In our case, we can see that the standard errors for all models are quite small, indicating that the mean estimates are likely quite precise.
Most Important Variables
The two most important variables for each model are consistently diabetes and arthritis across the board with regard to obesity. Specifically, people with diabetes and arthritis may be more likely to develop obesity.
Evaluating discrepancies: If a variable has a high importance in one model, but not in other models, it means that the other models are able to achieve comparable prediction accuracy without relying heavily on that particular variable. In this case, the variable BINGE is highly important in the lm model but not in LASSO, ridge, and enet models. This suggests that BINGE may be correlated with other predictors, and the regularization techniques in LASSO, ridge, and enet models are able to effectively deal with this multicollinearity by shrinking the coefficient of BINGE to zero or close to zero.
Therefore, it is likely that BINGE does not contribute much additional information in predicting obesity after accounting for the other predictors in the LASSO, ridge, and enet models. However, the lm model may be overfitting to BINGE due to the absence of regularization, which could explain why it is highly important in that model.
